EcoCommons_MarxanWeb_connection

Author

Zhao Xiang, EcoCommons

Introduction

Using the Species distribution modeling techniques provided by the EcoCommons Platform (www.ecocommons.org.au), we produced probability distribution maps for the three Queensland endangered species: koala, brush tailed rock-wallaby, and beach stone curlew.

Then we adjusted the probability distribution maps of these three species with the planning units shapefile prepared by the Marxan MaPP, and ran four planning scenarios with a target of expanding the coverage of protected areas in QLD to 30%.

EcoCommons Outputs

  1. Species records pulled from GBIF, ALA, EcoPlots, OBIS
  2. Species distribution modelling output: Species distribution Probability maps (This is the input tested in this project).

Marxan MaPP Inputs

  1. Shapefile of planning area and units.
  2. Shapefile of cost.
  3. Shapefile and csv of biodiversity features (Where EcoCommons can help!).

EcoCommons connects with Marxan Showcase:

  1. We get the QLD planning units from Marxan MaPP
#| echo: false

library(sf)
Linking to GEOS 3.12.2, GDAL 3.9.2, PROJ 9.4.1; sf_use_s2() is TRUE
library(terra)
terra 1.7.78
library(ggplot2)
library(ggspatial)

# We get the planning units shapefile via Marxan MaPP

QLD_Unit <- "/Users/zhaoxiang/Documents/tmp/ecocommons-marxan-integration-poc/ecocommons-marxan-integration-poc/QLD_Koala_Marxan/QLD_Unit/cost-surface-template_dbfd87df-99bc-4229-a607-04e632a67fda_20240828T015655.shp"

QLD_Unit  <- st_read(QLD_Unit)
Reading layer `cost-surface-template_dbfd87df-99bc-4229-a607-04e632a67fda_20240828T015655' from data source `/Users/zhaoxiang/Documents/tmp/ecocommons-marxan-integration-poc/ecocommons-marxan-integration-poc/QLD_Koala_Marxan/QLD_Unit/cost-surface-template_dbfd87df-99bc-4229-a607-04e632a67fda_20240828T015655.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 9778 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 137.904 ymin: -29.19145 xmax: 153.6074 ymax: -9.139518
Geodetic CRS:  WGS 84
# Get the number of rows
n_rows <- nrow(QLD_Unit)

# Plot the shapefile with no fill color and number of rows in the title
ggplot(data = QLD_Unit) +
  geom_sf(fill = NA, color = "gray") +
  theme_minimal() +
  ggtitle(paste("QLD Planning Unit", n_rows))

  1. I made a cost layer using the reciprocal of the distance to state-owned road as a surrogate of the cost.

    The assumption is: the closer to the state owned road, the more expensive to purchase the unit.

Reading layer `QLD_cost_road' from data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan/QLD_cost_road.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 9778 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 137.904 ymin: -29.19145 xmax: 153.6074 ymax: -9.139518
Geodetic CRS:  WGS 84

  1. Biodiversity features. I used EcoCommons to produce three species’ SDM to start with.

    Species 1: koala

    Species 2: brush tailed rock-wallaby

    Species 3: beach stone curlew

Loading required package: sp

  1. We need to turn these SDMs to binary results (shapefies).
#| echo: false

library(terra)
library(sf)


# Define a function to process each species
process_species <- function(raster_data, QLD_Unit, species_name, output_dir) {
  # Extract mean raster value for each polygon in QLD_Unit
  extracted_values <- extract(raster_data, vect(QLD_Unit), fun = mean, na.rm = TRUE)
  
  # Rename and update the feature column in QLD_Unit
  names(QLD_Unit)[names(QLD_Unit) == "cost"] <- "feature"
  QLD_Unit$feature <- extracted_values[,2]
  
  # Subset polygons where the feature (mean value) is >= 0.5
  QLD_species <- subset(QLD_Unit, feature >= 0.5)

  # Write the subset sf object to a shapefile
  shapefile_base <- file.path(output_dir, species_name)
  st_write(QLD_species, paste0(shapefile_base, ".shp"), delete_layer = TRUE)
  
  # List all the files that belong to the shapefile
  shapefile_files <- list.files(output_dir, pattern = paste0(species_name, "\\."), full.names = TRUE)
  shapefile_base_names <- basename(shapefile_files)

  # Zip the shapefile, only including the base filenames
  zipfile_path <- file.path(output_dir, paste0(species_name, ".zip"))
  old_wd <- setwd(output_dir)  # Change the working directory to the output directory
  zip(zipfile_path, shapefile_base_names)
  setwd(old_wd)  # Change back to the original working directory

  # Comment out the deletion step for now
  # file.remove(shapefile_files)
}

# Set the output directory
output_dir <- "/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan/"

# Apply the function to each species
process_species(koala, QLD_Unit, "koala", output_dir)
Deleting layer `koala' using driver `ESRI Shapefile'
Writing layer `koala' to data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan//koala.shp' using driver `ESRI Shapefile'
Writing 2905 features with 2 fields and geometry type Polygon.
process_species(beach_stone_curlew, QLD_Unit, "beach_stone_curlew", output_dir)
Deleting layer `beach_stone_curlew' using driver `ESRI Shapefile'
Writing layer `beach_stone_curlew' to data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan//beach_stone_curlew.shp' using driver `ESRI Shapefile'
Writing 1013 features with 2 fields and geometry type Polygon.
process_species(brushtailed_rockwallaby, QLD_Unit, "brushtailed_rockwallaby", output_dir)
Deleting layer `brushtailed_rockwallaby' using driver `ESRI Shapefile'
Writing layer `brushtailed_rockwallaby' to data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan//brushtailed_rockwallaby.shp' using driver `ESRI Shapefile'
Writing 782 features with 2 fields and geometry type Polygon.
### We can also turn them into a csv format

# Function to extract presence (1) and absence (0) from raster based on threshold (e.g., 0.5)
extract_presence_absence <- function(raster_data, unit) {
  extracted_values <- extract(raster_data, vect(unit), fun = mean, na.rm = TRUE)
  presence_absence <- ifelse(extracted_values[, 2] >= 0.5, 1, 0)
  return(presence_absence)
}

# Apply the function to each species raster
QLD_Unit$koala_presence <- extract_presence_absence(koala, QLD_Unit)
QLD_Unit$beach_stone_curlew_presence <- extract_presence_absence(beach_stone_curlew, QLD_Unit)
QLD_Unit$brushtailed_rockwallaby_presence <- extract_presence_absence(brushtailed_rockwallaby, QLD_Unit)

# Select relevant columns to create a presence-absence data frame
presence_absence_df <- data.frame(
  puid = QLD_Unit$puid,  # Assuming 'ID' is the unique identifier for the planning unit
  Koala = QLD_Unit$koala_presence,
  Beach_Stone_Curlew = QLD_Unit$beach_stone_curlew_presence,
  Brushtailed_Rockwallaby = QLD_Unit$brushtailed_rockwallaby_presence
)

# Write the presence-absence data to a CSV file
write.csv(presence_absence_df, file.path(output_dir, "presence_absence_species.csv"), row.names = FALSE)

# Check the CSV output
print(head(presence_absence_df))
  puid Koala Beach_Stone_Curlew Brushtailed_Rockwallaby
1    1     0                  0                       0
2    2     0                  0                       0
3    3     0                  0                       0
4    4     0                  0                       0
5    5     0                  0                       0
6    6     0                  0                       0
  1. Plot species SDM binary shapefile outputs for double check
library(sf)
library(ggplot2)


# Load the species shapefiles
koala_sf <- st_read(file.path(output_dir, "koala.shp"))
Reading layer `koala' from data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan/koala.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 2905 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 138.6178 ymin: -29.19145 xmax: 153.6074 ymax: -16.38993
Geodetic CRS:  WGS 84
beach_stone_curlew_sf <- st_read(file.path(output_dir, "beach_stone_curlew.shp"))
Reading layer `beach_stone_curlew' from data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan/beach_stone_curlew.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1013 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 138.0468 ymin: -28.45818 xmax: 153.6074 ymax: -9.139518
Geodetic CRS:  WGS 84
brushtailed_rockwallaby_sf <- st_read(file.path(output_dir, "brushtailed_rockwallaby.shp"))
Reading layer `brushtailed_rockwallaby' from data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan/brushtailed_rockwallaby.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 782 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 143.7571 ymin: -29.19145 xmax: 153.6074 ymax: -16.16685
Geodetic CRS:  WGS 84
# Add an identifying column to each species dataset for the legend
koala_sf$species <- "Koala"
beach_stone_curlew_sf$species <- "Beach Stone-curlew"
brushtailed_rockwallaby_sf$species <- "Brushtailed Rock-wallaby"

# Combine the species shapefiles into one dataset
combined_species_sf <- rbind(koala_sf, beach_stone_curlew_sf, brushtailed_rockwallaby_sf)

# Plot the unit (base map) first and overlay the species habitats without borders
combined_plot_with_unit <- ggplot() +
  geom_sf(data = QLD_Unit, fill = NA, color = "grey") +  # Base map (QLD Unit)
  geom_sf(data = combined_species_sf, aes(fill = species), color = NA, alpha = 0.6) +  # No borders for species
  scale_fill_manual(values = c("blue", "green", "coral")) +  # Colors for each species
  theme_minimal() +
  labs(title = "Species Habitats within QLD Unit", subtitle = "Koala, Beach Stone-curlew, Brushtailed Rock-wallaby") +
  theme(legend.title = element_blank())

# Display the plot
print(combined_plot_with_unit)

library(sf)
library(terra)

# Function to extract presence (1) and absence (0) from raster based on threshold (e.g., 0.5)
extract_presence_absence <- function(raster_data, unit) {
  extracted_values <- extract(raster_data, vect(unit), fun = mean, na.rm = TRUE)
  presence_absence <- ifelse(extracted_values[, 2] >= 0.5, 1, 0)
  return(presence_absence)
}

# Apply the function to each species raster and add to the QLD_Unit shapefile
QLD_Unit$koala_presence <- extract_presence_absence(koala, QLD_Unit)
QLD_Unit$beach_stone_curlew_presence <- extract_presence_absence(beach_stone_curlew, QLD_Unit)
QLD_Unit$brushtailed_rockwallaby_presence <- extract_presence_absence(brushtailed_rockwallaby, QLD_Unit)

# Save the updated QLD_Unit as a new shapefile
output_shapefile <- file.path(output_dir, "QLD_Unit_with_species_presence.shp")
st_write(QLD_Unit, output_shapefile, delete_layer = TRUE)
Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
Shapefile driver
Deleting layer `QLD_Unit_with_species_presence' using driver `ESRI Shapefile'
Writing layer `QLD_Unit_with_species_presence' to data source 
  `/Users/zhaoxiang/Documents/QCIF_doc/QLD_Koala_Marxan//QLD_Unit_with_species_presence.shp' using driver `ESRI Shapefile'
Writing 9778 features with 5 fields and geometry type Polygon.
# Verify the presence/absence columns
print(head(QLD_Unit))
Simple feature collection with 6 features and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 137.904 ymin: -17.06073 xmax: 138.475 ymax: -16.83687
Geodetic CRS:  WGS 84
  puid cost                       geometry koala_presence
1    1    1 POLYGON ((138.0468 -16.9487...              0
2    2    1 POLYGON ((137.904 -16.94876...              0
3    3    1 POLYGON ((138.0468 -17.0607...              0
4    4    1 POLYGON ((137.904 -17.06073...              0
5    5    1 POLYGON ((138.3323 -17.0607...              0
6    6    1 POLYGON ((138.1895 -17.0607...              0
  beach_stone_curlew_presence brushtailed_rockwallaby_presence
1                           0                                0
2                           0                                0
3                           0                                0
4                           0                                0
5                           0                                0
6                           0                                0

Marxan Four scenarios solutions:

Our SDMs input to Marxan MaPP:

EcoCommons SDMs output of three species on Marxan MaPP

Scenario 1: No SDMs, No Costs

No Costs, neither SDMs

Scenario 2: SDMS, No Costs

SDMs only

Scenario 3: Costs, No SDMs

Costs only

Scenario 4: SDMs + Costs

Costs and SDMs